##################
# Angel M. Villegas-Cruz
# This file replicates the main results, namely Figure 4, Figure 5, and Appendix Table 1.
# 08/18/2024
##################

library(ggplot2) 
library(stargazer)
library(httr)
library(tidyverse)
library(rworldmap)
library(dplyr)
library(sandwich)

options(scipen = 999)

mydata <- read_csv("mydata_final.csv") #upload data

#Analysis starts here

##################
#Appendix Table 1
##################

#Conducting logistic regression model 1
f1 <- glm(threat ~ ucdp_civil_war + muslim + elec_demo + 
            ln_gdppc +
            agree + ln_trade_share_i +
            ln_followers_count + ln_statuses_count + 
            ln_friends_count + xinjiang_num + factor(month),
          data = mydata, family = "binomial")
#Conducting logistic regression model 2
f2 <- glm(bene_rule ~ ucdp_civil_war + muslim  + elec_demo +
            ln_gdppc + agree + ln_trade_share_i +
            ln_followers_count + ln_statuses_count + 
            ln_friends_count + xinjiang_num + factor(month),
          data = mydata, family = "binomial")

#Calculating robust standard errors
cov1        <- sandwich::vcovHC(f1, type = "HC1")
robust_se1    <- sqrt(abs(diag(cov1)))
cov2       <- sandwich::vcovHC(f2, type = "HC1")
robust_se2    <- sqrt(abs(diag(cov2)))

#Creating Appendix Table 1
stargazer(f1,f2,type='text', out="appendix_table1.html",
          title="",
          dep.var.caption = "Repression image management strategies",
          dep.var.labels = c("Threat","Benevolent rule"),
          covariate.labels = c("Civil war",
                               "Muslim-majority",
                               "Electoral democracy",
                               "Log(GDP p.c.)",
                               "UN voting with China",
                               "Log(Trade share with China)",
                               "Log(N. of followers)",
                               "Log(N. of tweets)",
                               "Log(N. of friends)",
                               "N. of tweets about Xinjiang"),
          omit.stat = "f",
          se        = list(robust_se1,robust_se2), 
          add.lines = list(c("Month fixed-effects", "Yes","Yes")),
          column.sep.width = "-15pt",  font.size = "normalsize",digits=2)

#Note: Since they are too long, the main manuscript does not include all the month binary variables (i.e., month fixed effects).

#Open Appendix Table 1 in browser
BROWSE("appendix_table1.html")

##################
#Odds ratio plots (Figure 4 and Figure 5)
##################

# Create labels for the plots
boxLabels = c("Civil war",
              "Muslim-majority",
              "Electoral democracy",
              "Log(GDP p.c.)",
              "UN voting with China",
              "Log(Trade share with China)",
              "Log(N. of followers)",
              "Log(N. of tweets)",
              "Log(N. of friends)",
              "N. of tweets about Xinjiang")


##################
#Figure 4
##################

paste(f1$coefficients[2:11],sep="",collapse = ",") #Copy-paste coefficients below
paste(robust_se1[2:11],sep="",collapse = ",") #Copy-paste robust standard errors below
coeffici <- c(1.11824473751921,-0.0495804582974618,-0.412533122228517,0.0475980970222306,-2.08792765448439,0.0329581896697359,-0.149283677335546,0.0774003372589112,-0.140597379308877,0.00427922545658371)
std = c(0.395304880509999,0.24794868635087,0.270727966098656,0.0933349470525806,0.599151253367272,0.10369202193025,0.0561438639732976,0.134327418670017,0.0796509680181931,0.000967280572074619)

lci <- exp(coeffici - 1.96 * std) #calculate lower confidence interval
or <- exp(coeffici) #calculate odds ratio
uci <- exp(coeffici + 1.96 * std) #calculate upper confidence interval
lreg.or <- cbind(lci, or, uci)        
lreg.or

paste(lreg.or[,1],sep="",collapse = ",") #lci
paste(lreg.or[,2],sep="",collapse = ",") #or
paste(lreg.or[,3],sep="",collapse = ",") #uci

#Creating data frame for the new plot
df <- data.frame(
  yAxis = length(boxLabels):1,
  boxCILow  = lreg.or[,1],
  boxOdds = lreg.or[,2],
  boxCIHigh = lreg.or[,3]
)

# Creating Figure 4
p <- ggplot(df, aes(x = boxOdds, y = yAxis))
dev.new()
p + geom_vline(aes(xintercept = 1), size = .25, linetype = "dashed") +
  geom_errorbarh(aes(xmax = boxCIHigh, xmin = boxCILow), size = .5, height = .2, color = "gray50") +
  geom_point(size = 3.5, color = "blue") +
  theme_bw() +
  theme(panel.grid.minor = element_blank()) +
  theme(text = element_text(size = 15)) +
  scale_y_continuous(breaks = df$yAxis, labels = boxLabels) +
  scale_x_continuous(breaks = seq(0,7,1) ) +
  coord_trans(x = "log10") +
  ylab("") +
  xlab("Odds ratio (log scale)") +
  geom_text(label = "", y =1.1, x = 2.5, color="darkgray")+
  ggtitle("")

#Save as "Figure 4.pdf"

##################
#Figure 5
##################

paste(f2$coefficients[2:11],sep="",collapse = ",") #Copy-paste coefficients below
paste(robust_se2[2:11],sep="",collapse = ",") #Copy-paste robust standard errors below
coeffici <- c(1.96225370886596,1.04351783040401,0.763628402946896,-0.0737706477988885,-2.36895570515148,0.0326227182096438,-0.111709669410753,0.665774335292321,-0.472524177947535,0.00455768187885459)
std = c(0.273584453640462,0.195079145544509,0.224646130937335,0.0735866963753981,0.461342456076302,0.0781038547449544,0.0448755710812616,0.12188076679084,0.0679438983052492,0.000630520676231299)

lci <- exp(coeffici - 1.96 * std) #calculate lower confidence interval
or <- exp(coeffici) #calculate odds ratio
uci <- exp(coeffici + 1.96 * std) #calculate upper confidence interval
lreg.or <- cbind(lci, or, uci)        
lreg.or

paste(lreg.or[,1],sep="",collapse = ",") #lci
paste(lreg.or[,2],sep="",collapse = ",") #or
paste(lreg.or[,3],sep="",collapse = ",") #uci

#Creating data frame for the new plot
df <- data.frame(
  yAxis = length(boxLabels):1,
  boxCILow  = lreg.or[,1],
  boxOdds = lreg.or[,2],
  boxCIHigh = lreg.or[,3]
)

# Creating Figure 5
p <- ggplot(df, aes(x = boxOdds, y = yAxis))
dev.new()
p + geom_vline(aes(xintercept = 1), size = .25, linetype = "dashed") +
  geom_errorbarh(aes(xmax = boxCIHigh, xmin = boxCILow), size = .5, height = .2, color = "gray50") +
  geom_point(size = 3.5, color = "blue") +
  theme_bw() +
  theme(panel.grid.minor = element_blank()) +
  theme(text = element_text(size = 15)) +
  scale_y_continuous(breaks = df$yAxis, labels = boxLabels) +
  scale_x_continuous(breaks = seq(0,7,1) ) +
  coord_trans(x = "log10") +
  ylab("") +
  xlab("Odds ratio (log scale)") +
  geom_text(label = "", y =1.1, x = 2.5, color="darkgray")+
  ggtitle("")

#Save as "Figure 5.pdf"